공간정보분석
ASS-10.5
‘Der Kriging Krieg’
by: 서용훈
※ 본 문서는 PC 환경에 최적화되어있습니다.
1 미세먼지(\(PM_{10}\))의 분포
- 서울특별시를 대상으로
요즘 대두되고 있는 미세먼지를 주제로 삼아 Kriging 보간법을 사용하여 서울시의 미세먼지 분포도를 제작해보았습니다.
1.1 사용한 데이터셋
사용한 데이터셋으로는 2019년 1월부터 12월까지의 서울 도시 대기 확정자료와, 도시대기 측정망 관측소 위치자료, 마지막으로 서울시 셰이프 파일을 사용하였습니다.
1.2 자료 손질
1.2.1 사용한 패키지
자료손질시 사용하였던 패키지 목록입니다.
1.2.2 자료 불러오기
반복-조건문을 사용하여 12개월치 자료를 순차적으로 불러온 뒤, 한 데이터프레임으로 만들었습니다.
## data de air
fi <- paste0("./data/2019AQdata/", list.files(path = "./data/2019AQdata/", pattern = "확정자료엑셀다운*"))for (i in 1:12) {
assign(paste0("A", i, "raw"), read_xls(fi[i]))
if (i == 1) {
Araw <- plyr::rbind.fill(A1raw)
rm(A1raw) # 부산물 삭제
} else {
Araw <- plyr::rbind.fill(Araw, eval(as.symbol(paste0("A", i, "raw"))))
rm(list = as.character(paste0("A", i, "raw"))) #부산물 삭제
}
}
Araw 날짜 시도 측정소명 측정소코드 아황산가스 일산화탄소 오존
1 2019-01-01 01 서울 중구 중구 111121 .003 .7 .003
2 2019-01-01 02 서울 중구 중구 111121 .003 .8 .002
3 2019-01-01 03 서울 중구 중구 111121 .003 .9 .002
4 2019-01-01 04 서울 중구 중구 111121 .003 .8 .002
5 2019-01-01 05 서울 중구 중구 111121 .003 .8 .002
6 2019-01-01 06 서울 중구 중구 111121 .003 .8 .002
7 2019-01-01 07 서울 중구 중구 111121 .003 .7 .006
이산화질소 PM10 PM2.5
1 .054 39 24
2 .056 38 27
3 .057 42 28
4 .054 42 31
5 .048 49 33
6 .048 47 32
7 .04 48 30
[ reached 'max' / getOption("max.print") -- omitted 218993 rows ]
1.2.3 자료의 일 평균 구하기
시간대별로 수치가 기록되어있는 자료를 일평균을 내주는 작업입니다.
먼저 측정소 순으로 정렬한 뒤, YYYY-MM-DD HH로 이루어진 날짜 형식을 날짜와 시간으로 분리하여 일평균을 구할 수 있도록 날짜 열을 만들었습니다. 1
Araw$PM10 <- as.numeric(Araw$PM10)
order_Air <- Araw[c(order(Araw$측정소코드)), ]
Araw <- separate(data = Araw, col = 날짜, into = c("날짜", "시"), sep = " ")다음은 일 평균을 구하는 과정입니다.
for (j in 1:as.numeric(length(Araw$측정소코드)/24)) {
tmp <- vector(mode = "numeric", 24)
for (k in 1:24) {
if (is.null(Araw$PM10[(j - 1) * 24 + k]) == TRUE) {
davg[j] <- NULL
break
} else {
tmp[k] <- Araw$PM10[(j - 1) * 24 + k]
}
}
davg[j] <- mean(tmp)
}
rm(tmp)
dmean_air <- (Araw[c(seq(from = 1, to = length(Araw$PM10), by = 24)), ] %>% cbind(davg))[,
c(1, 3:5, 12)]
dmean_air 날짜 시도 측정소명 측정소코드 davg
1 2019-01-01 서울 중구 중구 111121 35.95833
25 2019-01-02 서울 중구 중구 111121 33.04167
49 2019-01-03 서울 중구 중구 111121 35.25000
73 2019-01-04 서울 중구 중구 111121 54.70833
97 2019-01-05 서울 중구 중구 111121 60.04167
121 2019-01-06 서울 중구 중구 111121 42.41667
145 2019-01-07 서울 중구 중구 111121 50.75000
169 2019-01-08 서울 중구 중구 111121 37.66667
193 2019-01-09 서울 중구 중구 111121 46.62500
217 2019-01-10 서울 중구 중구 111121 51.25000
241 2019-01-11 서울 중구 중구 111121 70.95833
265 2019-01-12 서울 중구 중구 111121 89.04167
289 2019-01-13 서울 중구 중구 111121 106.66667
313 2019-01-14 서울 중구 중구 111121 155.20833
337 2019-01-15 서울 중구 중구 111121 122.45833
[ reached 'max' / getOption("max.print") -- omitted 9110 rows ]
위처럼 반복문을 사용하여, 24개의 열씩 합하였습니다.
1.2.4 시간대 추출
해당 분석을 위하여 19년도 미세먼지를 검색하던 도중 눈에 띄는 뉴스를 하나 보게 되었습니다.
위 뉴스를 접하고 1월중 가장 미세먼지가 심한 날을 알아보았습니다.
for (j in 1:31) {
# 31일까지
if (isTRUE(dmean_air$davg[j] > max)) {
max <- dmean_air$davg[j]
date <- dmean_air$날짜[j]
} else {
max <- max
date <- date
}
}[1] 155.2083
[1] "2019-01-14"
신기하게도 뉴스가 방영된 날짜와 똑같음을 알 수 있었습니다.
해당 사례로 분석을 진행하기 위해 19년 1월 14일의 관측소 데이터만 추출하였습니다.
빈 데이터프레임을 만들고,
for (j in 1:length(dmean_air$davg)) {
if (dmean_air$날짜[j] == date) {
maxair[j, ] <- dmean_air[j, ]
}
}
maxpm <- filter(maxair, !is.na(davg))반복문을 사용하여 추출한 뒤, 결측치를 제거하는 방식으로 빈 열들을 제거하였습니다.
1.2.5 기존 자료와 주소 조인
해당 관측소의 위치를 알기 위하여 에어코리아에서 구득한 관측소 위치자료를 기존 데이터프레임에 조인하였습니다.
1.3 지오코딩
1.3.1 사용한 패키지
네이버 클라우드에서 제공하는 지오코딩 API 서비스를 사용하기 위한 API 관련 패키지입니다. 기존 구글로 지오코딩하는 방식은 지도의 최신화가 안되는 등의 이유로 국내 포털보다 정확도가 떨어져서 네이버 지오코딩 서비스를 사용하였습니다. 네이버 클라우드
1.3.2 지오코딩
우선, 발급받은 API ID와 비밀번호를 입력하였습니다. (API키는 개인 키라 가렸습니다.)
API 호출을 통해 지오코딩을 하는 과정입니다. 기존 API 호출로는 하나의 주소밖에 처리하지 못하므로, 이 또한 반복문을 사용하여 한번에 처리될 수 있게 하였습니다.
x <- vector(mode = "double", length(maxpm_refined$davg))
y <- vector(mode = "double", length(maxpm_refined$davg))
for (l in 1:length(maxpm_refined$davg)) {
apiResult <- httr::GET(url = "https://naveropenapi.apigw.ntruss.com/map-geocode/v2/geocode",
httr::add_headers(`X-NCP-APIGW-API-KEY-ID` = clientId, `X-NCP-APIGW-API-KEY` = clientSecret),
query = list(query = maxpm_refined$`측정소 주소`[l] #주소
) #리버스지오코딩시 쿼리 추가
)
if (apiResult$status_code == "200") {
result <- rawToChar(apiResult$content)
Encoding(result) <- "UTF-8"
list <- xmlToList(result)
u <- as.numeric(list$addresses$x)
v <- as.numeric(list$addresses$y)
x[l] <- u
y[l] <- v
} else if (apiResult$status_code == "400") {
print(paste("Bad Request Exception in", maxpm_refined$측정소명[l]))
} else if (apiResult$status_code == "500") {
print(paste("Unexpected Error in", maxpm_refined$측정소명[l]))
} else {
print("몰라요")
}
}
maxpm_refined$x <- x
maxpm_refined$y <- y1.4 Kriging
1.4.1 사용한 패키지
크리깅 보간 및 도식을 위해 사용한 패키지입니다.
1.4.2 서울시 행정경계
eTL에서 구득한 서울시 행정경계 파일을 마스킹을 통한 추출 및 CRS획득을 위해 불러왔습니다.
seoulpoly <- readOGR(dsn = "./data/shp", layer = "tl_emd._seoul_4326", encoding = "UTF-8",
stringsAsFactors = FALSE)OGR data source with driver: ESRI Shapefile
Source: "C:\Users\dydgn\OneDrive\바탕 화면\Analytical_Methods_for_Spatial_Information_Assignment\ASS-11\data\shp", layer: "tl_emd._seoul_4326"
with 467 features
It has 3 fields
1.4.3 교차검증을 위한 Samples
3개의 관측소 데이터를 무작위로 추출한 후, 분석하려는 데이터프레임에서 해당 자료를 제거하였습니다. 추출한 데이터는 추후 교차검증시 사용하였습니다.
1.4.4 공간자료로의 변환
지오코딩하여 좌표를 획득한 자료를 공간자료로 변환하였습니다.
1.4.5 경향면 구축
Kriging시 Detrending을 위한 경향면을 만들었습니다.
우선, 추후 레스터 자료의 출력을 위해 빈 격자형 자료를 만들었습니다.
library(gstat)
grd <- as.data.frame(sp::spsample(PMdatapts, "regular", n = 50000))
names(grd) <- c("x", "y")
coordinates(grd) <- c("x", "y")
gridded(grd) <- TRUE
fullgrid(grd) <- TRUE
proj4string(grd) <- proj4string(PMdatapts)1st Order
1차식으로 만든 경향면입니다.
# Define the 1st order polynomial equation
f.1 <- as.formula(davg ~ x + y)
# Run the regression model
lm.1 <- lm(f.1, data = PMdatapts)
# Use the regression model output to interpolate the surface
dat.1st <- SpatialGridDataFrame(grd, data.frame(var1.pred = predict(lm.1, newdata = grd)))
r <- raster(dat.1st)
plot(r)2nd Order
2차식으로 만든 경향면입니다.
# Define the 2nd order polynomial equation
f.2 <- as.formula(davg ~ 1 + x + y + I(x * x) + I(y * y) + I(x * y))
# Run the regression model
lm.2 <- lm(f.2, data = PMdatapts)
# Use the regression model output to interpolate the surface
dat.2nd <- SpatialGridDataFrame(grd, data.frame(var1.pred = predict(lm.2, newdata = grd)))
r <- raster(dat.2nd)
r.m <- mask(r, seoulpoly)
tm_shape(r.m) + tm_raster(n = 10, palette = "-RdBu", auto.palette.mapping = FALSE,
title = "Predicted PM10") + tm_shape(PMdatapts) + tm_dots(size = 0.2) + tm_legend(legend.outside = TRUE)위와 같이 1차와 2차식간 차이가 명확하여, Kriging시 2차 경향면을 변수로 삼았습니다.
위에서 나온 잔차를 속성프레임에 집어넣어 추후 분석에 사용하기 쉽도록 하였습니다.
1.4.6 Variogram
Kriging 보간을 위해 필수적인 Variogram을 제작하였습니다.
# library(automap) fit.vgm = autofitVariogram(prodp ~1, tzprice1_v, model =
# 'Exp') plot(pr.v, fit.vgm$var_model)변수 바이오그램
fvgm <- fit.variogram(emvgm, vgm(psill = 100, "Hol", range = 15, nugget = 0.5), fit.method = 6)
plot(emvgm, fvgm, xlim = c(0, 30))잔차 바이오그램
잔차 바이오그램을 만들던 도중,
위와같이 맨 왼쪽에 이상치 하나가 있어서, 보다 나은 fit을 위해 이상치를 제거하였습니다.
fvgm.res <- fit.variogram(emvgm.res, vgm(psill = 100, "Exp", range = 5, nugget = 0)) #resid vario
plot(emvgm.res, fvgm.res, xlim = c(0, 30))1.4.7 Ordinary Kriging
이번 분석에서는 Ordinary, Simple, Universal 총 3가지 크리깅 기법을 사용하였습니다.
위 기법들 중 Ordinary Kriging을 실시하는 과정입니다.
잔차
[using ordinary kriging]
PM10 값
[using ordinary kriging]
CV
다음은 교차검증과정입니다. 추출한 표본의 수가 적어 의미가 적으므로, 참고만 부탁드립니다.
cv V2
1 164.5076 157.9583
2 168.7822 171.0833
3 160.8667 150.6667
se <- vector(mode = "numeric", 3)
for (m in 1:3) {
se[m] <- (PMdataCVpts$davg[m] - cv[m])^2
}
rmse <- sum(se)/3 %>% sqrt()
rmse[1] 87.88937
Var Surface
크리깅 연산에서 중요한 부산물인 분산도를 아래와 같이 제작하였습니다.
1.4.8 Simple Kriging
특정 위치에서 지역화돠된 변수에 대한 기대값은 전체지역에서의 평균값과 같다고 가정하는 단순 크리깅을 시도하였습니다. 보통 정규 크리깅을 많이 사용하는데, 그 이유는 단순 크리깅에서 지역화된 변수에 대한 가정이 현실적이지 않기 때문입니다.
PM10 값
[using simple kriging]
CV
cv <- extract(RSK, PMdataCVpts)
cvlm <- lm(cv ~ PMdataCVpts$davg)
plot(cv ~ davg, data = PMdataCVpts) cv V2
1 161.6311 157.9583
2 168.8942 171.0833
3 161.9564 150.6667
se <- vector(mode = "numeric", 3)
for (m in 1:3) {
se[m] <- (PMdataCVpts$davg[m] - cv[m])^2
}
rmse <- sum(se)/3 %>% sqrt()
rmse[1] 84.14248
Var Surface
1.4.9 Universal Kriging
크리깅 연산시 경향을 고려한 크리깅입니다.
PM10 값
[using universal kriging]
CV
cv <- extract(RUK, PMdataCVpts)
cvlm <- lm(cv ~ PMdataCVpts$davg)
plot(cv ~ davg, data = PMdataCVpts)
Call:
lm(formula = cv ~ PMdataCVpts$davg)
Residuals:
1 2 3
0.5902 -0.2108 -0.3794
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 140.76599 8.01776 17.56 0.0362 *
PMdataCVpts$davg 0.13421 0.05007 2.68 0.2273
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.7326 on 1 degrees of freedom
Multiple R-squared: 0.8778, Adjusted R-squared: 0.7556
F-statistic: 7.185 on 1 and 1 DF, p-value: 0.2273
se <- vector(mode = "numeric", 3)
for (m in 1:3) {
se[m] <- (PMdataCVpts$davg[m] - cv[m])^2
}
rmse <- sum(se)/3 %>% sqrt()
rmse[1] 102.32
Var Surface
2 마치며
저번 강의에서 배운 보간법인 크리깅을 미세먼지에 적용켜보았으나, 생각보다 그리 잘 들어맞지 않았던 것 같습니다. 미세먼지는 공기라는 유체에 떠다녀서 다른 변수들 또한 고려를 해야하는데 이 부분이 부족했던 것 같습니다. 이러한 부분을 보완하기 위하여 공동크리깅(Co-Kriging)을 사용하여야 할 것 같다고 사료됩니다.